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SECTION I 
INTRODUCTION 


This study considers the development of a digital 4D automatic 
control law to capture a steep glides I ope, flare and land a CTOL a i r- 
craft under low visibility and adverse wind conditions using the 
Microwave Landing System (MLS) under development by the FAA Cl], C2H, 

The study of curved 4D flight paths leading to a steep final approach 
under low visibility (up to CAT MIC) and adverse wind conditions is 
part of the Terminal Configured Vehicle (TCV) program sponsored jointly 
by NASA and FAA [3]. The goals of the TCV program include the reduction 
of aircraft noise in airport communities, the reduction of the adverse 
effects of weather conditions on the efficiency of aircraft operations 
in air terminals, the reduction of fuel consumption, and the efficient 
use of airspace in congested terminal areas through the use of the 
Microwave Landing System and the development of technology. 

The major effect of the use of steep gli deslopes in the final 
landing is in the area of noise reduction. in comparison to the 
currently used 2.5° to 3° l LS glideslopes, a 5° to 6° glideslope reduces 
the noise level on the ground due to differences in the altitude pro- 
files and thrust levels. At equal distances from the runway, the 
altitude of an aircraft on a 6° glideslope is almost twice the altitude 
of an aircraft on a. 3° glideslope. Thus, if both aircraft were to 
generate the same amount of noise, the noise level perceived on the 
ground would be reduced due to the 2:1 ratio in the altitudes. A 
further reduction in noise is due to the fact that the aircraft 
following a 6° glideslope requires a lower thrust setting than the one 
following a 3° glideslope. For a Boeing 737 the reduction in the thrust 
levels is of the order of 2:1. Thus, for a steep glideslope less 
engine noise is generated. 

The reduction in thrust level associated with steep glideslopes 
also reduces the fuel consumption during the final approach. The 
ability to fly glideslope of varying angles may also provide a method 
to avoid the vortex generated by large aircraft, by allowing smaller 
aircraft to fly different glideslopes to reduce the likelihood of such 



encounters. In general, the ability to fly steep glides I opes provides 
a versatility that can be useful in efficient use of airspace in the 
terminal area. 

The guidance information necessary to fly steep glideslopes under 
low visibility conditions is provided by the Microwave Landing System 
(MLS). The MLS is a ground-based guidance system providing position 
information to aircraft which are inside its volumetric coverage at 
discrete intervals of time. Thus, aircraft have accurate position in- 
formation to fly glideslopes of varying angles, and follow curved 4D 
paths to landing within the MLS coverage under low visibility conditions. 
The MLS consists of a DME providing range information, an azimuth 
antenna co- located with the DME providing the aircraft T s azimuth angle 
relative to the runway, and an elevation antenna located at the glide- 
path intercept point (GPIP) but offset to the side of the runway pro- 
viding the aircraft with its elevation angle relative to the local 
horizon. The azimuth antenna provides coverage up to + 60°, the eleva- 
tion coverage is up to 20°. A second elevation antenna located further 
down the runway is also under consideration for flare guidance. Thus, 
the MLS provides position information with a high accuracy in a volumet- 
ric coverage, with low sensitivity to adverse weather conditions. 

In this study, a longitudinal guidance and control law using MLS 
data and body-mounted acce I erometers is developed for the final approach 
phases from glideslope capture to touchdown. 

In Section II, the aircraft equations of motion and the wind model 
are described. A mathematical model describing the perturbations of 
the aircraft from a 6° glideslope is obtained for the longitudinal 
axis; the first order effects of lags in the effects of winds on the 
motion of the aircraft and thrust build-up are included in this model. 
Using the Dryden spectra, a dynamic model of gusts is developed; these 
steady winds and wind shear models are combined to obtain a total wind 
model. These mathematical models are used in simulating the aircraft's 
motion under various wind conditions. In Section III, the desired 
flight path is described and formulated in a form suitable to the 
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development of the digital control law. In Section IV, the development 
of the digital control law is described. The control law consists of 
a filter estimating error parameters and wind velocities which are used 
as feedback to compute control commands. Section V describes the 
results obtained from a digital simulation of the aircraft using the 
guidance and control law developed. 
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II. AIRCRAFT DYNAMICS AND WIND MODELING 


The aircraft equations of motion are used in simulating aircraft T s 
flight and in the development and analysis of the control law. Similarly, 
the wind model is used in generating various turbulence levels and the 
development of the filtering equations. To obtain a suitable mathematical 
model of the aircraft dynamics and the effects of various wind conditions 
on the motion of the aircraft, the general equations of motion can be 
linearized about a steep glideslope at an airspeed of 120 knots with 
fully extended flaps and landing gear down. Thus, a linear mathematical 
model describing the perturbat ions of the aircraft motion from the 
steady flight condition of glideslope tracking is obtained. This model 
can then be used to simulate the aircraft T s motion and for analysis. 

In the general perturbation equations (e.q. C4], pp. 2.33), the 
aircraft's lateral and vertical motions are coupled. This coupling is 
mainly due to a non-zero roll angle. As the phases of flight considered 
are glideslope capture, glideslope tracking and flare (as shown in 
Figure I), the roll angle is small, and can be expressed as a pertur- 
bation from a steady value of zero. In this case, the lateral and 
vertical motions are decoupled. Thus, the longitudinal equations of 
motion can be used to simulate the aircraft's motion in the vertical 
plane. The longitudinal equations can then be expressed in the stability 
axes as 


mu = - mg cosy Q 9 + f + fj 

a x x 


( la) 


m(w - U q) = - mg sin y 0+f A + f_ 
o 0 A z ^z 


(lb) 


I q = m A + m-r ( I c) 

yy M A T 


where 


U q - steady inertial 


speed 


i n the x d i rection 
s 
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"Y 


gl ides I ope angle 

perturbation in pitch angle 
pitch rate 

perturbation in net aerodynamic force along the x direction 


perturbation in net aerodynamic force along the z direction 


perturbation in thrust along the x direction 


perturbation in thrust along the z^ direction 


perturbation in pitching moment due to aerodynamic forces 
perturbation in pitching moment due to thrust 


moment of inertia about the y axis 

1 s 


a i rcraf t mass 


u - perturbation in inertial speed along the x^ direction 
w - inertial speed along the z s direction 

Mainly, three sets of coordinate axes are used in this study: the 

earth fixed axes, the body axes, and the stability axes. The earth-fixed 

coordinate frame ( x ^, y^, z ) has its origin fixed at the glidepath 

intercept point (GPIP) on the runway. The x^ axis is along the runway 

centerline, the direction in which the aircraft lands being chosen 

positive velocity along x . The z^ axis is the local vertical positive 

downwards; y is perpendicular to both x and z , with its positive end 
e e e 

directed so as to make the coordinate frame right-handed. The earth is 
assumed to be stationary with respect to inertial space; so that the 
earth-fixed axes form an inertial frame. 


The body axes (x , y , z ) and the stability axes (x , y , z ) are 

D D D S S S 

fixed to the body of the aircraft; i.e. they are body-fixed axes. The 
origin of both axes is fixed at the aircraft center of mass. The axis 
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is along the fuselage reference line of the aircraft, positive towards 
the nose, the axis is positive towards the tip of the right wing; the 
axis is perpend i cu I ar to both x fe and and is positive downwards (when 
the aircraft pitch angle is zero). The stability axes (x g , y g , z^) are 
obtained from the body axes by a rotation of ct Q , the steady state angle of 
attack, about the y^ axis. Three sets of coordinate frames are shown in 
Figure 2. 

In equation (I), the term f^ represents the total algebraic 

x 

change in the value of the aerodynamic force along the x^ axis due to 

changes from steady state values in the values of the aerodynamic and 

control variables; the terms f^ , f^ , f , m^ and m_p are defined 

z z x 

similarly as the changes in the appropriate forces or moments from their 
steady values on the glideslope. These terms can be expressed in terms 
of the aircraft stability derivatives, the moments of inertia and the 
perturbat i ons in the aerodynamic and control variables. Substituting 
these expressions for the force and moments into equation (I) the air- 
craft equations of motion can be expressed as given below. 


mu = 


= - mg cosy 0 + q Q S{(-C + 2C^ + C^. + 2C^_ ) u 1 

° u u . T o xu ? xo 


+ (C Lo ” So 5 - " C DS Se ~ C DSS 6s + C TXST 6 T} 
— e 


(2a) 


m(w -u q ) = - mg s i ny 0 + q S {-(C. . + 2C. ) u ? 

o H y 1 o ^o Lu ? Lo — 


(C 


La 


+ C n ) 
Do 




_q} 


(2b) 


I q = q Sc {(C., , + 2C.. ) u ? + (C M + C.._ ) a + C • a 
yy M ^o Mu T Mo — MTa_ — Ma, — 


+ C.. q + C.. x 6e + C MX <5 C + C..« T <5T}. 
Mq -1 M6e M6s S M6T 


(2c) 


where 


u f 


u + u 


(3a) 
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a = a + a , (3b) 

— w 

q = q + q w , (3c) 

u , a , and q are the components due to wind, q is the steady 
WWW o 

value of the dynamic pressure, S is the effective wing area, 5e, 5s 
and 5T are the perturbat i ons of elevator, stabilizer and thrust, 
respectively. 

Equation (2) describes the linear and angular velocities of the 
aircraft. The position of the aircraft can be obtained by integrating 
the inertial velocity components over time. Thus, the 4D error in 
position (i.e., the actual position relative to the desired position at 
a given time) can be linearized and expressed in the following form: 

• • 

x - x 

= - siny 0 + cosy u' + siny a (4a) 

Uo o o .o 

• • 

z - z 

— Uj— = - cosy^ 0 - siny^ u T + cosy^ x (4b) 

where x (t) and z (t) are the x and z coordinates of the desired 
o o e e 

aircraft position on the glideslope. The desired glideslope position 
in this case was taken to be the position of the aircraft going down 
a 6° glideslope at 120 knots inertial speed. 

The- effects of the servo responses of the actuators were modeled by 
first and second order systems. Thus, the thrust and stabilizer models 
are 


6T = - .5 ST + .298 5th, , 

(5a) 

6th = U3 , 

(5b) 

5s = us ; 

(5c) 
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where 5T is the thrust perturbation in thousands of pounds, 6th is the 
throttle position perturbation in degrees and 6s is the stabilizer 
position perturbation in radians. As the lags in elevator action are 
small, the elevator time constants were neglected. The aircraft 
equations of motion (2), the position equations (4) and the actuator 
equations (5) can be combined and after some manipulation expressed in 
state variable form as shown in equation (6). 

X - X z - z 

x' = (9 u' a q — 2. — ST s+h 5s), (6a) 

^ Uo Uo 

u ? = (6e 6s 5th), w T = (u ! a q ) (6b) 

WWW 

x = Ax + Bu + Dw , (6c) 

where A, B and D are matrices of appropriate size, corresponding to the 
coefficients in equations (2), (4) and (5). These matrices are given 
in terms of the aircraft stability derivatives in the Appendix. For 
a more detailed derivation of the models given here see [5H. Thus, 
the aircraft equations of motion and the 4D error equations for sink 
rate and ground speed have been modeled in state variable form as shown 
i n equation (6) . 

The wind model contains the components of steady wind velocities, 

turbulence and shear winds in the longitudinal axis. The turbulence 

model uses the Dryden spectra C6H for the various components varying 

with altitude. The turbulence model has three components: in the 

x, direction, a in the z, direction and q which models the effects 
b g b H g 

of turbulence on the pitch rate of the aircraft. These components are 
modeled using the following spectra. 

2 a 2 L 

S (Q) = — , (7a) 

u I + (L fl) 2 

u 

o 2 l | + 3(L fl) 2 

s (n) = — — - — , (7bj 

“ Va 2 [| + (L fl)2]2 
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S (Q) = 

q 


fi 2 Va 2 


+ (i5. n) 2 « 

7T 


s (n) 


(7c) 


where b is the wing span, and L w are the scales of turbulence, Va is 
the airspeed, and Q is the spatial frequency related to the temporal 
frequency w by 


G = w/Va . 


(8) 


The u T component is independent of a and q , however, a and q are 
9 9 9 9 9 

correlated with their cross spectral density being 


S (w) = 
9 « 




, . 4b 

+ J — rr~ w 
J irVa 


S (w) 
a 


(9) 


The above spectra can be factored using spectral f actor i zat ion 
methods to obtain a linear system driven by white noise which generates 
an output having the above spectral characteristics C7H, [8]. Thus, 
the following transfer functions are obtained to generate u T ^, a g an( ^ 9g- 


G (s) = 
u 




G (s) = 
a 


+ /3 

Va 


+ + (^) 2 s 2 

Va Va 


( 10a) 


( 10b) 


G (s) = 

q 


, 4b 
uVa S 


( 10c) 


where a is the input to the system G (s) to obtain q with the 
9 9 9 

specified spectrum and cross-spectra I density. Figure 3 shows a block 

diagram of the system generating the turbulence components. 
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The steady and shear wind in the longitudinal direction was 
modeled by 


u* = u 1 j u? , = W 3 , 

s sg sh 5 


( I la) 


a = W 4 
s 4 


(Mb) 


Thus, to simulate a specified shear profile for u T s , with 
appropriate initial conditions; e.g., to obtain a linear profile 
u ! s changing at a rate of u f s ^ q , the initial condition for is 

set to u f shQ and W 3 is set equal to zero; alternately, an impulse in 
W 3 will also achieve the same profile. 

The transfer functions obtained for the wind model can equivalently 
be expressed as differential equations in state variable form as shown 
in (12); the matrices A^, and are given in the Appendix. 

W = A w W + B w 5 , w = C w W , (12) 

where w is given by equation ( 6 b). 
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III. FLIGHT PATH MODELING 


In this section, the desired flight path for the aircraft will be 

modeled as a vector differential equation. In Section IV, the solution 

to the optimal control problem of following a flight path specified in 

state variable form will be presented. The flight path of interest 

here consists of a steep glideslope and flare. First consider the 

flight path describing the glideslope, which will be called the reference 

flight path as the equations of motion obtained in the preceeding 

section are linearized about this flight condition. Let be a vector 

where represents the pitch angle, Z^ the speed in the direction, 

Z the speed in the z direction, Z the pitch rate, Z and Z the 
r 3 s r 4 r 5 r 6 

x q and z q components of the aircraft position respect i ve I y , Z^ the 

thrust force, Z the throttle position, and Z the stabilizer 
r 8 r 9 

position. Then, the reference glideslope flight path can be described 
by the following set of differential equations. 


Z = Z 

r i r 4 


Z =0 
r 2 


z = ° 

r 3 


Z =0 
r~4 


(13) 


Z =- U siny Z + cosy Z + siny Z + £ 
r 5 o o ri 'o r 2 o r 3 £5 

Z rg =-U o cosy 0 Z ri - sin Yo Z rz + cosy o ^ 

Z =0 

r? 


Z r = 0 
r 8 


Z =0 
r 9 


where c =U 0sinY ,5 = U © c° s Y^ 

^5 o o o’ -^Tg o O o 


(14) 
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is the steady state pitch angle, and the initial conditions are 
given by 


Z 

r l 



Z 

r 2 




(15a) 


^ r 6 




z 

r 9 


S 

o 


(15b) 


Thus, the glidepath condition can be described by the equation 


A Z 
— r 


+ G. 


(16) 


where the matrix A_ corresponds to the coefficients in equation (13). 

The desired flight path can also be described by an equation similar to 
equation (16). 


Zi - Z 4 + 

Z2 = £2 

Z 3 *53 (,7) 

Z 4 = £4 

Zs = - U siny Zi + cosy Z2 + siny Z 3 + £5 
3 o 0 1 0 o — 

Z 6 = - u o cos Y 0 z i - siny o Z 2 + cosy o Z 3 + £ 6 



To define the desired flight path completely, £. must be specified. 
Thus, let and be the desired position coordinates in the earth- 
fixed coordinate frame, then setting 

* • • • __ 

£2 ~ Xj cosCZj - a Q ) - z d sin(Zi - a Q ) - Z3Z4 

£3 = *x d sin(Zj - a Q ) + *z’ d cosCZj - a ) + Z2Z4 (18) 

£5 x d [I - cos(Z x - Z ri )] + z d sin(Zi - Z^) + U o 'sir>Y 0 Zj 

£6 =-x^ sin(Zi - Z pi ) + z d Cl - cos(Zi - Z^)3 + U q cosy q 

would result in 

Z 5 = x d > Z 6 = Z d ' 

while Z2 and Z 3 would correspond to the desired velocity components in the 
stability axes. Note that £1 , £4, £7, £9 and £9 should be chosen 
appropriately, so as to aid in the performance of following the desired 
flight path; these are specified with due cons i derat ion to the aero- 
dynamic characteri st ics of the aircraft. Thus, the desired flight 
path can also be expressed in the form of equation (16). 

Z = A Z + £ (19) 

Subtracting the reference flight path Z^ from the desired flight 
path Z, we obtain the deviations of the desired flight path from the 
reference glidepath. Further, normalizing the velocity components 
Z2 - Z r2* Zs “ , and the position components Z5 - Z^, Zg - Z^ by 

the nominal speed of the aircraft, we form a normalized vector, say z, 
describing the desired flight path relative to the reference flight 
path. Using (13) - (19), z can be expressed as 

z = A z z + c , (20) 


where 
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A = 
z 


0 

0 

0 

0 

-sinY c 

-cosy c 

0 

0 

0 


0 

0 

0 

0 

COSY 


0 

0 

0 

0 

sirw 


-siny cosy 


0 

0 

0 


0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


(21 ) 


• • • • 

*d z h 

? 2 = 0 - cos(Y o + z x > - g- sin(Y 0 + z^ - Z3Z4 
o o 


• • 

X 


£3 = si n(Y Q + z i> + [P c°s(Y q + z\ ) + (I + Z 2 ) Z4 
o o 


x d Z d 

£5 = g— ( I - cosz^ + g—sinzx+sinY zi 

o o ° 


( 22 ) 


X d Z d 

£g = g — s i nz x + g — ( I - coszx ) + cosy Zx 
o o J 


Since the aircraft dynamics equations (6) are expressed as 
deviations from the reference glidepath, z describes the desired trajec- 
tory values for the vector x in equation (6). Thus, in the following we 
shall cal! z the desired flight path. Now, to define the desired 
flight path z, it is necessary to specify £, hence x^, x^, z d and z^. 

When the aircraft is following the glideslope, the desired 
deviation from the reference flight path is zero. Thus, an initial 
condition of zero for all the components of z and setting £ equal to 
zero in equation (20) generates the glideslope as the desired flight 
path. During flare, the desired flight path differs from the glideslope; 
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I 


During this portion of the flight z wi I I thus be different than zero. 
The flare flight path was determined by a constant deceleration in the 


vertical and longitudinal directions. Thus, the desired flare path 

starts with a specified sink rate z, (i.e. the glideslope sink rate), 

* y • • 

and a specified ground speed x^, and using constant decelerations 

and decelerates to specified speeds z^ and at a specified point 

z df and x df and then maintains its speed till touchdown. With the above 

constraints, it is possible to determine the unspecified variables to 


obtain the following relations 



x - x 
df dg 


(23a) 


(23b) 


* x +d 

Z df 


(23c) 



(23d) 


(23e) 
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Thus, specifying the desired velocity at the start of flare (i.e. 

• • • 
x, > z, ), the desired velocity at the end of flare (i.e. - z, x ), and 
dg dg 1 df 

the desired touchdown point (i.e. x^) the remaining variables can be 
determined using equation (23). The instaneous velocity components 
during flare can be computed using 

x (t) = x - X* (t - t ) , t > t (24a) 

d dg d g - g 

z,(t) = z, + Y,(t = t ) , t > t (24b) 

d dg d g - g 

where t is the time at which flare is initiated. 

g 

Using equation (23) and (24), x^, z^, x" d and z’ d can be computed, 
then these values are used to compute instantaneous values for c; in 
equation (22) which generates the desired flight path z given by 
equation (20). The desired pitch during flare was generated by a 
constant pitch rate until the desired pitch angle reached a specified 
value 0^. The thrust, throttle and stabilizer positions were set to 
their values on the glideslope; however, the throttle and stabilizer 
components were not weighed in the cost function, thus allowing 
deviations from their desired values without any associated penalty. 

This allowed the use of throttle and stabilizer to be determined by the 
dynamics of the aircraft required to follow the desired flight path. 

Thus, the approach taken here allows the desired flight path to 
be specified in a simple manner by choosing the velocity and accele- 
ration associated with a given path. The desired flight path was de- 
scribed by a set of differential equations in state variable form as 
given in equation (20). Thus, the values generated by equation (20) 
for the state vector z are the values that the aircraft dynamics state 
vector should ideally have. Hence, the difference between the actual 
aircraft state and the desired state, i.e. x - z, describes the error or 
the deviations from the desired flight path, the goal of controlling the 
aircraft being to keep these deviations as small as possible. 
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IV. DEVELOPMENT OF FILTER AND CONTROL LAW 


In the previous sections, the aircraft’s longitudinal motion was 
modeled by equation (6) including the response to control surface and 
thrust changes under various wind conditions. The wind condi-tions 
were modeled by steady winds, random gusts and shear in equation (12), 

In the last section, the desired flight path during glideslope tracking 
and flare was expressed by equation (20). In this section, these 
models are combined to develop a control law to follow the desired 
flight path; i.e. to capture a steep glideslope, follow the glideslope 
until flare altitude and follow the flare path to touchdown under 
various wind conditions including gusts and shear winds, utilizing MLS 
data and body-mounted acce I erometer data. 

The control law is developed by imbedding the problem of following 
a specified flight path into a quadratic regulator with disturbance 
problem [9], DO]. As the separation principle holds for this problem, 
the control law consists of a filter in cascade with a feedback control 
law. 

A. Development of the Control Law 

The desired flight path consists of capturing the glideslope, 
following the glideslope til! flare altitude, then following the flare 
path to touchdown. This flight path is described by equations (20) 

(23) and (24) in state variable form. The difference between x and z 
is the error in the actual flight path. Combining equation (6) and (20), 

x - z = Ax + Bu + DW - A^z - £ : (25) 

rearranging terms, 

e = Ae + Bu + Dw + (A - A^) z - £ , (26a) 

e = x - z . (26b) 
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The last three terms in the right-hand-side of equation (26a) can 
be combined as follows. 


w r d = (w T z f s T ) , 

D , = (D (A - A,) -I) 

d z 

e = Ae + Bu + 0 ,w , 
d d 


(27a) 

(27b) 

(27c) 


From equation (27c), we can see that the problem of following a 
desired trajectory can also be interpreted as a regulator problem with 
disturbances by considering the error in the actual trajectory; i.e., 
the difference between the actual and desired tra jector i es . To use the 
results for this problem, the disturbances must be expressed in state 
variable form and the cost function expressed in terms of the error, e. 
Thus, combining the wind model with the model for the desired flight 
path, we get 


d = d + n , 


w d = c d d » 


where 


d f = (w T z T £ r ) , 


A J = 


0 

0 


0 

A 


0 

0 

0 


C, = 


0 

I 


0 

0 

I 


(28a) 

(28b) 


(29a) 


(29b) 


(29c) 
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n 


(29d) 


B 


w 


e 



Note that if we let be a gaussian white noise process, then £ 
would be a Brownian motion process. Thus, z(t) would be statistically 
described by a spectrum shaped by equation (20). Using this approach, 
the desired flight path is assumed to be a sample function (or a 
realization) of a random process z(t) with a specified power spectral 
density. The control law would be independent of the specific sample 
function, and be applicable to a class of flight paths. An important 
advantage of this method over the usual quadratic tracking problem 
(Cl G pp. 219-227) is the fact that the desired flight can be specified 
on line, i.e. the desired flight path can be altered without having to 
make any modifications to the control law. A second advantage of this 
method is that a set of constant gains can be obtained by solving for 
the steady state solution, which would not be possible in a finite time 
prob I em. 

Since the objective is to follow the desired flight path z, 
minimizing a cost function which is quadratic in the error x - z is 
appropriate; hence, 

, f f 

J = J E / C(X - 2 ) ' Q (x - z) + u' R ull dt . (30) 

o 

Substituting equation (26b) for the error term, we obtain 
. + f 

J = j E / He* Q e + u* R uD dt , (31 ) 

o 

where E is the statistical expectation operator. Now, since a digital 
control law is desired, the control variables can vary only at distinct 
invertals of time say at sampling intervals of T seconds, thus, 

u (t) = u R , kT < t '< KT + T - (32) 
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Thus, equations (27c), (28), (31) and (32) together with linear 
measurements pose a quadratic regulator problem with random disturbances 
which has been treated in [5H, C9]. With equation (32), equation (27c) 
can be viewed as a sampled-data system and can be expressed in discrete 
form [_5j 

e k+ i = + e k + r u k *■ r d d k + n k ’ (33) 

where t is the transition matrix correspond i ng to A, r is the system 
control gain, is the system disturbance gain and is a white noise 
sequence. The disturbance equation (28) can also be expressed in 
discrete form 


k + 


►d d k 


+ n, 


(34) 


Where is the transition matrix correspond i ng to A^, and n^ is 
a white noise sequence. The cost function J can also be equivalently 
expressed in discrete form using equations (27) and (28), E5U, to 


obta i n 

J i = 2 E J e 'k^ e k + u 'k R u k + 
k=o 


T 

0 = f f (s) Q (j> (s) ds , 
o 

. T 

R = j ri (s) Q r (s) ds + RT , 
o 

. T 

M = / <f> T (s) Q r (s) ds , 
o 

T 

N = / 1 (s) 0 r d (s) ds , 

o 

T 

S = J rf d (s) Q r (s) ds 
o 


2 (e f . Nd. + e\Mu. + d T .sH 
k k k k k 

(35) 

(36a) 

(36b) 

(36c) 

(36d) 

(36e) 
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Note that the continuous and discrete cost functions J and Jx 
differ by a constant which does not depend on the control hence for 

minimization purposes the two cost functions are equivalent, i.e., the 
same control minimizes both J and Jx« Thus, the continuous optimal 
control problem can be reduced to a discrete one; i.e. minimize Jx 
subject to the constraints of equations (33) and (34) while the control 
u k may only use past measurements. The solution to the continuous 
problem posed by equations (27c), (28), (31) and (32) as well as the 
discrete problem posed by equations (33), (34) and (35) is given by the 
following control law C5D. 

% - - H > k ; k - % k ■ (37 > 

Where e^ and d^ are the least mean-square estimates of and d^ 
given past measurements respect i ve I y , and the gains and are 
defined by the equations 

H lk = R k 2 Gl k ' H2 k = R51 Gz k ’ 

Gi k = r Pi k + m , G 2k = r D k + s , 

D k* p >k rd + %*<> - \ p >k r 

while Pj and P 2 k are given by the nonlinear difference equations: 

p l k _ ( = *’ p l k + - R- 1 G lk + Q , Pi N = 0 , (39a) 

k 

p 2 k _| = C* - r R-l G lk ]' D k + N , P 2n = N . (39b) 

For the finite time problem where t f and N are finite, the control 

gains vary with time; however, under very loose restrictions C5H, the 
gains converge to constant values. These constant gains are used in the 
control system. It should also be noted that the steady state gains 


(38a) 

(38b) 

(38c) 
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result in a stable closed loop system. To see the effect of the various 
components, we can write the control as 


u 

k 


CH e «k + H » "k + H z 2 k + H c ^ 


(40) 


From equation (40), we see that the control law uses the estimate 
of current flight path error e^ to reduce this error to zero; it further 
uses wind estimates to prevent expected disturbances due to winds 
including shear winds; the last two terms provide lead control to follow 
the desired flight path. Note that contains velocity terms as well 
as position. While contains acceleration information. Thus, the 
control law consists of a filter which receives measurements obtained by 
the MLS and on-board sensors followed by the feedback term given in 
equation (40). A functional block diagram of the system is shown in 
Figure 4. 

B. Development of Fi Iter Equations 

The measurements of aircraft position are obtained from the 
Microwave Landing System (MLS) [2H, which provides azimuth, elevation 
and range to the aircraft at discrete intervals of time; the other 
aircraft variables measured by on-board sensors include pitch, pitch 
rate, barometric altitude, barometric sink rate, calibrated airspeed, 
and body mounted acce I erometer data. In terms of the aircraft 
variables defined so far in this study, the measurements were expressed 
as fo I lows. 


Y l = 0 O + x i + V X 


(41a) 


Y 2 = x 4 + V 2 


(41b) 


2 2 |/2 

Y 3 = [Oj x 5 + x - x) +(u x 6 + z - z)H +V 3 


Yu = tan -1 ( 


IL x 6 + z - 2 




U x 5 + z 
u o 3 o 


(41c) 

(4ld) 
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(4 1 e) 


= - (U o X 6 


+ Z + Vc) 


y 6 = - (U o + z Q ) V 6 

I + x 2 

Y 7 = U [x 2 tan x 3 + x 3 

O o 

C03 z X 3 


( 4 1 f ) 

x 4 (I + x 2 + COS X ! -COSX 3 )] + V 7 

(4lg) 


I + u' 


Yo = U. (- 


o cos a 


-> Vc 


Y 9 = U q Lx 2 + (I + x 2 ) x 4 tan x 3 + x 3 sin x. 


(4th) 

x 4 sin x x ] *+ V 9 (4 I i ) 


where x. is the component of the aircraft state vector x defined by 

equation (6a), x and z are the current 4D coordinates of the desired 

glideslope, z is the vertical coordinate of the elevation antenna 
© 

phase center, x and z are the coordinates of the azimuth antenna 
a a 

phase center in the earth-fixed coordinate system and V. is the noise 
t h 1 

introduced by the i T sensor. The earth-fixed coordinate system is 

referenced to the point where a perpend icu ) ar from the elevation antenna 

phase center crosses the runway centerline. The expressions for the 

accelerations along z s and x g , Y 7 and Yg respect i ve I y , were obtained by 

differentiating the correspond i ng velocities in the earth-fixed 

coordinates and transforming them back to the stability axes. The 

slant range, Y 3 , and the elevation, Y 4 , were obtained for a planar MLS. 

The values for the standard deviations of the sensor noises were chosen 

to reflect current i nstrumentat ion standards £\2j, II 1 3U ; the values 

are given in Table I. The measurements given by equation (41) are 

then processed by a non-linear transformation given by equation (42) to 

obtain flight path error signals which are fed into a filter. 


(42a) 

(42b) 


(]“ (r e cos y 4 + x q ) - z 5 = e 5 + v 3 


Yl = Y i “ 0 O " Z 1 = ® x + vj 
y 2 = Y 2 - z 4 = e 4 + v 2 


V3 = ~ 


(42c) 
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(42d) 


Vh = - Jj~ ( p e sin Y 4 " z e + Zq ) - z 6 = e 6 + v 4 
o 

Y5 = -Q— <z Q - Y 5> - z 6 = e 6 + v 5 

O 

I T f 

Y6 = — (siny o - Y 6 ) - E 6 Az = E 6 Ae + v 6 

o 

y 7 = - E3 (Bu + Az) = E3 (Ae + DC w W) + v 7 

u o 

y 8 = ■gS. - (I + Z 2 ) = e 2 - Ej C w W + v 8 
o 

y 9 = -g 2 - - E 2 Az = E 2 (Ae + DC w W) + v g 
o 


where 


r 

e 

n = 


= CY 2 - L 
3 


2 

e 


h 


tan" 1 


ea 

x 

a 


2L 

e 


Y 3 cos (Y 4 + 


x 

a 

cos n ' 


( 42e) 
( 42f ) 
(42g) 
(42h) 
(421) 

(43a) 

(43b) 


h is the height of the elevation antenna above the azimuth antenna, 
ea 

r is the ranqe of the aircraft from the elevation antenna, E. is the 
e J 

column of the identity matrix, and v. is a pseudo noise representing 

J 

the additive noise in the measurement y . . Thus, we see from equation 
(42) that y can be expressed in terms of the error vector e and the wind 
components W. 


y k = y (kT) = Ci e k + C 2 W R + V k 


(44) 


where Ci and C 2 correspond to the coefficients of e and W in the right- 
hand-side terms in equation (42). Thus, y k is a measurement vector 
which can be expressed as a linear combination of e and W. Now, 
rewriting (26a) and (12), 
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(45a) 


e = Ae + DC W + Bu + (A - A,)z - £ 
w z 

w = A w w + B w 5 , (45b) 

it is seen that the flight path error, e, and the wind components, W, 
can be estimated using a Kalman filter. Since the desired flight path 
is known z and £ are not estimated but used as known parameters. 

Treating the problem as shown in C5H, and using the constant Kalman-Bucy 
ga i ns, we obta i n 


+ r u . + r 

k z 


z k - s k 


- C 2 W, ) (46a) 

^ k 


= W k + Fz (y k ■ C l e k ' ° 2 V 


(46b) 


where e^ and are filtered estimates of e^ and respectively given 
past values of y^. Thus, the filter provides the estimates of flight 
path error, e^, and wind estimates, W^, which along with the desired 
flight path parameters, and are fed into the control gain 
equation (40) to obtain the control commands, u^. A block diagram of 
the filter is shown in Figure 5. A block diagram of the control law 
comprising the filter and feedback is shown in Figure 6. 
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V. RESULTS 


To evaluate the performance and capabilities of the control law 
and filter developed in previous sections for a CTOL aircraft, a 
digital simulation of a Boeing 737 aircraft was developed. The flight 
path considered consisted of capturing a steep glideslope (6°) from 
level flight at 120 knots, following the glideslope to flare altitude 
and performing flare until touchdown with a nominal touchdown sink rate 
of 2 ft/sec. Thus, initially the aircraft is in a trim condition for level 
flight at 120 knots at an altitude of 3,153 ft., aligned with the runway, 
at a range of 30,000 ft. (I knot = 0.5144 m/sec; I ft. = 0.3048 m). 

The simulation includes the aircraft dynamics, wind gusts, steady 
winds, wind shear, filter, control law, and measurement generation. 

The aircraft dynamics were simulated using the perturbation model 
developed in Section II using the aerodynamic character i sties of a 
Boeing 737 aircraft. The wind conditions were generated using a random 
gust model with the Dryden spectra (including the effects of gusts on 
aircraft attitude rates), steady winds, and wind shear using the wind 
model developed in Section II. From the simulated aircraft position, 
velocity, attitude and attitude rates the on-board measurements were 
obtained using the measurement model given by equation (41) with 
additive or multiplicative sensor errors according to the sensor. The 
measurements used consist of MLS data at a rate of 10 samples per 
second, airspeed, barometric altitude and sink rate, attitude, attitude 
rates and body-mounted accelerometer data. These measurements are 
processed by a nonlinear preprocessor and then a linear filter with 
constant coefficients as shown in equations (42) and (46) respectively. 

The filter outputs the .estimates of flight path errors and wind 
velocities, which are input to the control gain box which generates 
elevator position, stabilizer rate, and throttle rate commands which 
are used in controlling the aircraft and to simulate the next point. 

Thus, the total closed-loop system is simulated. 

The simulation developed has the capability of generating various 
wind conditions, by adjusting the gust levels, the steady wind velocity 
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and the rate of change of wind velocities for wind shear generation. 
Sensor noise levels are also input parameters which can be altered. 

Logic is included in the filter to accommodate MLS dropouts accompanied 
by MLS flags. The start and stop times of such dropouts are also 
input parameters. During a MLS dropout, the filter goes into a mode 
where position is updated without using MLS data while other data is 
still used. Furthermore, initial flight path errors at the start of 
glides I ope capture and the sampling rate are also input parameters 
that can "be changed to simulate the closed loop response to different 
cond i t ions . 

Simulation runs of flights starting at glideslope capture and 
ending at touchdown were made under various conditions; these are shown 
in Figures 7-13. To compare the effects of various wind conditions, 
sensor errors, and MLS dropouts, a run with the ideal conditions of no 
measurement noise and no winds was made; this is shown in Figure 7. The 
aircraft motion for this case is smooth; at glideslope capture the a i r- 
craft pitches down and increases its sink rate to acquire the 6° 
glideslope. The flight path errors settle to zero following the glide- 
slope, which the control surfaces and thrust settle to their steady 
values for a 6° glideslope. The flight path errors settle to zero 
following the glideslope, while the control surfaces and thrust settle 
to their steady values for a 6° glideslope putting the aircraft into a 
trim condition for the glideslope. At flare the aircraft pitches up, 
reduces its sink rate and reduces its forward speed slightly as com- 
manded by the desired flight path. To avoid a high sink rate (as 
required by a steep glideslope) at low altitudes, the flare altitude 
was set at 200 ft. although flare initiation at lower altitudes provides 
acceptable flight path performance and touchdown conditions; it is 
expected that this provides a sufficient period of time and enough 
altitude for the pilot to react and take over should a failure to 
engage the flare rqode occur; however, a higher flare altitude can 
easily be obtained by changing the desired flight path parameters. It 
should also be noted that engaging flare at 200 ft. allows a more 
gradual and smooth flare maneuver for passenger comfort than a lower 
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flare altitude would allow. Note that during flare thrust is gradually 
increased, so that a higher engine pressure ratio is obtained at lower 
a f 1 1 tudes. 

The flight path followed under no wind and no measurement errors 
shown in Figure 7, is considered to be the nominal flight path as this 
is the response of the closed loop system under no d i sturbances . Thus, 
the nominal touchdown point is 1646 ft. past the runway threshold at a 
sink rate of 1.93 ft/sec. Thus, touchdown under various wind conditions 
is compared to these nominal values. However, the 4D horizontal and 
vertical errors plotted in Figures 7 to 13, are deviations from the 
desired flight path rather than the nominal flight path shown in 
Figure 7. The case with no winds but with measurement noise present 
is shown in Figure 8, where the presence of noise is seen to cause 
small deviations in the pitch axis. Figure 9 shows the closed loop 
response to a steady tail wind of 8 ft/sec. and average gust conditions 
( r =2 ft/sec., £ = I ft/sec.). It is seen that the wind velocity 

estimates follow the actual wind velocities; particu lar I y , the total 
wind velocity has little error. With the steady wind present, the 
control system retrims the aircraft to follow the desired flight path 
as can be seen from the pitch attitude. 

To see the closed loop response to a wind shear condition, under 
average gust conditions and errors in the initial estimates, a wind 
shear condition where a tail wind of 8 ft/sec. (about 5 knots) 
changes linearly into a 16 ft/sec. (about 10 knots) head wind from 
500 ft. altitude to 200 ft. was generated; this condition is shown in 
Figure 10. It is important to note that, for steep glideslopes, the 
effect of wind shear on the flight path accuracy is larger than for 
shallow glideslopes. This is due to the difference in sink rate at 
constant speed; a steep glideslope results in a higher sink rate, so 
that the aircraft goes from 500 ft. altitude to 200 ft. in a shorter 
period of time, and the same change of wind velocity occurs in a 
shorter period of time. Thus, as far as controlling the aircraft in 
wind shear is concerned, a 5 knot per hundred feet wind shear on a 
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6° glideslope has similar effects as an approximately 10 knot per 
hundred feet wind shear on a 3° glideslope. As can be seen from 
Figure 10, the wf nd' estimates follow the changing wind, and the control 
law using this information keeps the flight path errors low; in fact 
it is difficult to notice that a- shear condition is present by looking 
only at the 4D horizontal and vertical errors; however, control action 
and changes in the attitude point out the shear condition. It is also 
interesting to note that at the start of the wind shear, the filter 
attributes most of the change in wind velocity to. fast changing gusts, 
but as the wind velocity pers i sts in the change, the wind shear 
estimate follows the change. The total wind in the forward direction, 
however, is estimated more accurately on an instantaneous basis than 
its components, and has no noticeable lag; thus, this information can 
be displayed to the pilot in real time to provide quick warning of 
chang i ng wi nds. 

To see the effects of an MLS dropout condition, a 10 second 
dropout condition characterized by bad flags output by the on-board MLS 
receiver was simulated; this case is shown in Figure II. Thus, from 
30 sec. to 40 sec. the filter operates in a mode which does not use 
MLS data but updates position estimates using the other sensor data 
such as body-mounted acce I erometer and air data. The wind conditions 
correspond to. a tail wind of 8 ft/sec. and average gust levels. It is 
seen from the plots that the filter can update position without 
significant flight path error for a 10 second period during an MLS 
fade. Figures 12 and 13 show the cases where the gust levels are 
high = 4 ft/sec., 5^ = 2 ft/sec.) with a tail wind, with and without 
a wind shear condition. Table 2 shows the touchdown condition for the 
various runs. 
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VI. SUMMARY AND CONCLUSIONS 


In the proceeding sections, a digital automatic longitudinal 
control law for steep glideslope capture, flare to touchdown was 
developed using guidance data from the Microwave Landing System. The 
system uses MLS data, body-mounted accel erometer data, air data, 
attitude and attitude rate data to compute the control surface commands 
for automatic landing. The control law is developed to accommodate 
turbulence and wind shear whose effect is magnified for steep glides! opes. 
The filter used filters out the noise in the received data and computes 
estimates of the wind velocities which are used to generate control 
commands under the estimated wind conditions to follow the desired 
f I i ght path. 

A digital simulation of the aircraft and the control law was 
developed for use in a partial evaluation of the system developed. 

The simulation includes the longitudinal aircraft dynamics for a 6° 
glideslope, wind conditions in terms of gusts, steady winds and wind 
shear, sensor outputs with noise, filter and controller. As a partial 
evaluation of the performance of the system a limited number of runs to 
simulate various conditions were made. The flight path tracking per- 
formance evaluated by flight path errors show only small deviations on 
the glideslope and quick response in capturing the steep glideslope; 
the touchdown conditions show no excessive sink rates correspond i ng to 
"hard" landings, and no excessive touchdown dispersion. The wind 
estimates computed by the filter follow the wind velocities and changes 
in the wind velocities; wind shear conditions are detected through 
large and lasting changes in the wind velocity estimates. Thus, the 
simulation results indicate that the techniques presented can be used 
to perform automatic steep approaches with MLS guidance. For a more 
complete evaluation of the system developed, use of non linear aircraft 
dynamics coupled with a lateral control law for crab and decrab 
maneuvers would be required. 
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Figure l 1C Closed Loop Response Simulation: Meas. Noise Present, 
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Table 2 


Fig. No . 
and Wind 
Conditions 



No Winds 


Touchdown pt. 
( f t . ) 


1569 


Tail Wind 
Average Gust 


1503 


10 

Wind Shear 
Average Gust 


1607 


I I 

Wind Shear 
Average Gust 
(MLS Dropout 
on G1 i des I ope 


1607 


12 

Tail Wind 
High Gust 


13 

Wind Shear 
High Gust 


ft. 


0.3048 m 


Cond i t ions 


Sink Rate 
( ft/sec. ) 

Ground Speed 
( ft/sec. ) 

1 .93 

199.2 

2.63 

199.2 

2.28 

199.3 

1 .27 

199.2 

1 .27 

199.2 

2.83 

199.5 

1 .60 

199.2 
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APPENDIX 


The aircraft equations of motion which are used in the simulation 
and for the development of the filter and control law were developed 
in Section II. The final form of the continuous equations as given in 
equation (6c) is repeated here for convenience. 

x = Ax + Bu + Dw 

The form of the matrices A, B and D are also given in Section II. 
In this appendix, we shall give the expressions for the non-zero 
elements of these matrices in terms of the aircraft stability deriva- 
tives, the stability derivatives are assumed to be in the stability 
axis. 


a l 


mg cos 0 

q S 
M o 


a 2 


mg sin 0 ^ 
q S 


mU 

a 3 = — — , o 4 = 1/(03 + C, .•) , 

q S 


“5 = V^o SE ' L = i p u o 2 

where S is the wing area, c the mean aerodynamics chord and p is the 
air density. Using these variables, the matrix elements are given below. 


a 14 " 1 

°1 | 

a 21 = — , a 2 2 = — (-C n - 2 CL + C_ + 2C_ .) 

03 a3 Du Do Txu Txo 

-CL + C. C-r C nx 

Dot Lo Tx _ D6s 

a 2 3 ^ ’ 327 " — , 928 TT- 
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a 3 1 

= a 2 «4 , 

Q) 

CO 

ro 

II 

1 

«4 (C LU 

+ 2C Lo ) ' 

a 3 3 “ 

- *4 (C Lq 

+ C Do> • 

a 34 

= 

(a 3 

-•w 

» a 37 

" “4 C j Zr a 39 - 

" “ 4 C L6s 


a 41 

= a 4 

C • 

ma 

a 2 /a 5 , 

a 42 = 

— (C + 
as mu 

2 C 

mo 

+ C ,r,a 3 3 2 5 

y 

a 4 3 

= \__ 
«5 

(C + c _ 

ma mTa 

+ C • a 
ma 

33^ ' a 44 

- _L 

“5 ‘ 

(C + c . 

mq ma 

a 34> > 



c T + c • 

mT ma 

a 37 

c . + 

mos 

C • a 4 q 
ma 



a 47 

a 5 


' 349 “ a 5 




a 77 

= - .5 


, a 78 = .298; 




b 2i 

C D6e 

b 31 

= " 014 C L6e ' 

k C m6e 

+ C • b 3 x 
ma 


«3 

°41 

<*5 

y 

d 33 

= a 4 (C L . - 

C. ) 
Lq 

C + 

ril 3 - mq . . 

(d 33 - 1) 

C • 

ma 


> a 43 

*5 




Thus, given the stability derivatives of the aircraft, the A, B 
and D matrices can be computed. 
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